%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------  Simulate Equilibrium  --------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Setup workspace
close all
clear 

machine = '/Users/giulia/';

set(0,'defaultTextInterpreter','latex')

M = 15;

g_euler = 0.57721;

scale = 10^5;

%% Load matching function estimates

load(path_mf)
options = optimset('Algorithm','interior-point','display','on','MaxIter',90000,'MaxFunEvals',90000,'TolFun',10^-9,'TolX',10^-9,'TolCon',10^-9);

cobb_par = zeros(15,2);

for i =1:15
    col = freights(:,i)~=0;
    s = ships(col,i);
    f = freights(col,i);
    m = matches(col,i);
    lb = [0,0];
    ub = [Inf,0.95];
    theta0 = [rand(1),rand(1)];
    [cobb_par(i,:),~,~] = fmincon(@(theta)cobb(theta,f,s,m),theta0,[],[],[],[],lb,ub,[],options);
end

s_values = [1,linspace(0.95,1.95,10)];
v_values = linspace(0.75,5,11);

clearvars -except s_values v_values scale machine cobb_par

%% Load cost estimates

load(path_ships);

load(path_exports);

%% Frictionless allocation

Niter = 100000;
eps = 0.01;
diff = 10;
Winit = 0;
W0 = sum(w_eff(1:8));    
diff0 = Inf;

%
% Initialize the starting point
%

s_i = s_eff;
q_ij = m_eff.*Gdest_eff;
q_i = nansum(q_ij,2);
radius = 0.0;
delta_0 = (1 + unifrnd(-radius,radius,15,1)).*(q_i./s_i).*nansum(Gdest_eff.*meanprices_eff,2);
Gdest_0 = q_ij./nansum(q_ij,2);
lambda_0 = nansum(q_ij,2)./s_i;

%
% Setup convergence parameters
%

alpha = 10^-5;
coeff_q = 10^-6;
coeff_s = 10^-6;
iter = 0;
W_save = 0;
Niter = 10000000;

%
% We iterate over the matrix q_ij of carries traveling loaded on each route
% (ij) and over the number s_i of carriers waiting empty
%

while (iter<Niter)&&(diff>exp(-3))
    
    %
    % given a guess for q_ij and s_i, we compute the destination matrix G, 
    % matching rate lambda, and matches q_i
    %
    
    Gdest_0 = q_ij./nansum(q_ij,2);
    lambda_0 = nansum(q_ij,2)./s_i;
    q_i = nansum(q_ij,2);
    
    %
    % compute the average payoff delta_1 and value functions that implement
    % q_ij and s_i
    %

    [U_k,W_k,J_k,delta_1,inner_iter] = inner_loop(s_i,q_i,Gdest_0,lambda_0,delta_0,alpha,ksi,cs,cu,sigma,g_euler,S_obs,M,beta);
    
    % if the algorithm to compute delta_1 doesn't converge, choose a
    % smaller smoothing parameter

    if inner_iter == 10000
        alpha = alpha./10;
        [U_k,W_k,J_k,delta_1,inner_iter] = inner_loop(s_i,q_i,Gdest_0,lambda_0,delta_0,alpha,ksi,cs,cu,sigma,g_euler,S_obs,M,beta);
        [alpha, inner_iter] %#ok<NOPTS>
    end
    

    %
    % compute the gradient of the monopolist's profits with resepct to q_ij 
    % and s_i
    %
    
    [dpi_dq] = Dwelfare_Dq(q_ij,beta,delta,s_i,cobb_par,E,M,r,K_obs,c_inv,W_k,U_k,J_k,sigma,sf);
    [dpi_ds,~,~] = Dwelfare_Ds(q_ij,beta,delta,s_i,cobb_par,E,M,r,K_obs,c_inv,cu,sigma,W_k,U_k,J_k,sf);

    %
    % updte the guesses
    %

    coeff_mat = coeff.*ones(size(q_ij));
    coeff_mat(q_ij<0.01) = coeff_z;
    q0_ij = q_ij +coeff_mat.*dpi_dq;
    s0_i = s_i +coeff_s.*dpi_ds;
    
    
    iter = iter+1;
    q_ij = q0_ij;
    s_i = s0_i;
    delta_0 = delta_1;
end

q_star = nansum(q_ij,2);
Gdest_star = q_ij./nansum(q_ij,2);
ccp_star = ccp_equil(W_k,M);
s_star = s_i;
f_star = inverse_m(q_star,s_star,cobb_par,M,sf);
lambda_f_star = nansum(q_ij,2)./f_star;
etilde_ij = entrants(q_ij,s_i,cobb_par,M,delta,sf);
Gunc_star = zeros(size(Gunc_ineff));
Gunc_star(:,2:16) = etilde_ij./E;
Gunc_star(:,1) = 1-nansum(Gunc_star,2);
U_star = U_k;
W_star = W_k;
meanprices_star = prices(q_ij,beta,delta,s_star,cobb_par,E,M,r,K_obs,c_inv,sf);
Jf_star = zeros(M,M);
for i=1:M
    
    aa = (1 - delta.*beta.*(1-lambda_f_star(i)));
    
    Jf_star(i,:) = lambda_f_star(i).*(r(i,:)-meanprices_star(i,:))-c_inv(i,:);
    
    Jf_star(i,:) = Jf_star(i,:)./aa;
    
end
w_star = Welfare(Gdest_star,Gunc_star,ccp_star,s_star,f_star,q_star,cu./sigma,cs./sigma,c_inv,K_obs,ksi,r,E,sigma);
